//      mp_main.cc
//      matching puirsuit algorithms
//      
//      Copyright 2012 tqlong <tqlong@espada>
//      
//      This program is free software; you can redistribute it and/or modify
//      it under the terms of the GNU General Public License as published by
//      the Free Software Foundation; either version 2 of the License, or
//      (at your option) any later version.
//      
//      This program is distributed in the hope that it will be useful,
//      but WITHOUT ANY WARRANTY; without even the implied warranty of
//      MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
//      GNU General Public License for more details.
//      
//      You should have received a copy of the GNU General Public License
//      along with this program; if not, write to the Free Software
//      Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
//      MA 02110-1301, USA.

/* 
./main.template --help
 */

#include <iostream>
#include <boost/program_options.hpp>
#include <boost/date_time/posix_time/posix_time.hpp> 

#include <armadillo>

#include "helper/stl_helper.h"
#include "helper/ublas_helper.h"
#include "helper/ml_data.h"
#include "helper/optimization.h"

namespace po = boost::program_options;
po::options_description desc("Allowed options");
po::variables_map vm;

void process_options(int argc, char** argv);
void usage(int exitcode = 0);

void read_data();
void least_square();

#define NOW boost::posix_time::microsec_clock::local_time()
#define TIME_FROM_START (NOW - __start_time)
boost::posix_time::ptime __start_time = NOW;

int main(int argc, char** argv)
{
  //MSG("Matching pursuit ... ");

  process_options(argc, argv);

  if (vm.count("help"))
    usage();
  else if (vm.count("read"))
    read_data();
  else if (vm.count("ls"))
    least_square();
  else
    usage();

  MSG("Duration = " << TIME_FROM_START);
}

void process_options(int argc, char** argv)
{
  desc.add_options()
    ("help", "produce this usage")
    ("input", po::value<std::string>(), "input file")
    ("output", po::value<std::string>(), "output file")
    ("debug", po::value<std::string>(), "debug file")
    ("verbose", "print debug information")

    // GENERAL param

    // OPTIMIZATION param
    ("optim_iter", po::value<int>(), "optim: maximum number of iterations")
    ("optim_tol", po::value<double>(), "optim: tolerance")
    ("optim_mem", po::value<int>(), "optim: L-BFGS memory")
    ("optim_eta", po::value<double>(), "optim: fixed step size")

    // MODULE param

    //JOB type
    ("read", "read data file\n--input [--output]")
    ("ls", "Least Squares\n--input [--output]")
    ("test", "test UBLAS")
    ;
  po::store(po::parse_command_line(argc, argv, desc), vm);
  po::notify(vm);
}

void usage(int exitcode)
{
  std::cout << desc << std::endl;
  exit(exitcode);
}

void read_data()
{  
  RequiredParameter(vm, "read")("input");
  std::string inputFile = vm["input"].as<std::string>();
  typedef DataTable<sparse_vec_t, sparse_matrix_t> Data;
  //typedef DataTable<vec_t, matrix_t> Data;

  Data D;
  DataReader::read_libsvm(inputFile, D, true); // add a bias term
  MSG("sz(D)=" << sz(D) << " dim(D)=" << D.dim());
}

// square error function
struct LeastSquareFunction {
  struct var_t {
    vec_t x_;
    double b_;

    var_t(const vec_t& x, double b) : x_(x), b_(b) {}
    var_t add(const var_t& v, double s) const { return var_t(x_+s*v.x_, b_+s*v.b_); }
    var_t reverse() const { return var_t(-x_, -b_); }
    double dot(const var_t& v) const { return ublas::inner_prod(x_, v.x_) + b_*v.b_; }
    double norm() { double norm = ublas::norm_2(x_); return sqrt(sqr(norm)+sqr(b_)); }

    var_t add(const var_t& v) { return var_t(x_+v.x_, b_+v.b_); }
    var_t& addin(const var_t& v) { x_ += v.x_; b_ += v.b_; return *this; }
    var_t& addin(const var_t& v, double s) { x_ += s*v.x_; b_ += s*v.b_; return *this; }
  };
  typedef vec_t helper_t;

  matrix_t& A_;
  vec_t& y_;
  var_t init_x0_;

  LeastSquareFunction (matrix_t& A, vec_t& y) 
  : A_(A), y_(y), init_x0_(zero_vec_t(A_.size2()), 0) {}
  
  // compute residual
  void helper(const var_t& v, helper_t& r) {
    r = ublas::prod(A_, v.x_) + scalar_vec_t(A_.size1(), v.b_) - y_;
  }
  
  double operator() (const var_t& v, const helper_t& r) {
    double norm_r = ublas::norm_2(r);
    return 0.5*sqr(norm_r);
  }

  void grad(const var_t& v, const helper_t& r, var_t& grad) {
    grad.x_ = ublas::prod(ublas::trans(A_), r);
    grad.b_ = ublas::sum(r);
  }
};

void least_square()
{
  RequiredParameter(vm, "ls")("input")("optim_iter");
  std::string inputFile = vm["input"].as<std::string>();
  typedef DataTable<vec_t, matrix_t> Data;

  Data D;
  DataReader::read_libsvm(inputFile, D, false); // add a bias term
  MSG("sz(D)=" << sz(D) << " dim(D)=" << D.dim()); 

  LeastSquareFunction ls(D.data(), D.label());
  GradientDescent<LeastSquareFunction> gd(ls, VM(vm, "optim_iter", int));
  
  LeastSquareFunction::var_t x0(ls.init_x0_);
  gd.search(x0);

  MSG("weight=" << gd.x_.x_ << " bias=" << gd.x_.b_);
}
